Correlation energies beyond the random-phase approximation: ISTLS applied to 

spherical atoms and ions 
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The inhomogeneous Singwi, Tosi, Land and Sjolander (ISTLS) correlation energy functional of 
Dobson, Wang and Gould [PRB 66 081108(R) (2002)] has proved to be excellent at predicting 
correlation energies in semi-homogeneous systems, showing promise as a robust 'next step' fifth-rung 
functional by using dynamic correlation to go beyond the limitations of the direct random-phase 
approximation (dRPA), but with similar numerical scaling with system size. In this work we test the 
functional on spherically symmetric, neutral and charged atomic systems and find it gives excellent 
results (within 2mHa/e~ except Be) for the absolute correlation energies of the neutral atoms tested, 
and good results for the ions (within 4mHa/e~ except B + ). In all cases it performs better than the 
dRPA. When combined with the previous successes, these new results point to the ISTLS functional 
being a prime contender for high-accuracy, benchmark DFT correlation energy calculations. 
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Since their development, density- functional theory [1, 
2] (DFT) methods have vastly increased the range of 
quantum mechanical problems that can be studied. This 
wide range comes through the use of approximation 
to the exchange-correlation (xc) physics necessarily in- 
troduced to make Kohn-Sham (KS) theory possible. 
The most common approximations such as the LDA[2], 
GGA[3] and hybrid schemes [4] perform generally well, 
but usually give very poor results for electron correla- 
tion alone. In particular they give completely incor- 
rect physics for van der Waals (vdW) dispersion physics, 
which governs the weak bonds between widely sepa- 
rated systems. This physics can be important in sys- 
tems where there is a realm of near-zero density between 
sub-systems, such as in stretched molecules or lattices. 
The vdW physics is reintroduced in the popular vdW- 
DF[5] group of functionals, however such methods fail to 
reproduce the correct exponent [6] for vdW power laws 
U = —C p D~ p in zero-gap systems with at least one long 
and one short dimension, such as thin slab geometries 
and nano-wires[7]. 

An alternative approach to total energy calculations 
is to: i) solve for a groundstate under a given scheme 
(e.g. LDA) to evaluate V (r) and, via the KS Hamil- 
tonian h = — ^V 2 + V KS (r), to evaluate the orbitals and 
KS energies through hipi = tiipi and the density through 
n ( r ) = J2i /i|V'i( r )| 2 where fi is the occupation number 
of orbital i; ii) recalculate the energy using the so-called 
exact exchange (EXX) functional for exchange and a dif- 
ferent functional for correlation. Here we use the Hartree 
and exchange pair density ri 2 Hx(f , r') = n(r)n(r') — 
I J2i fii } i( r )' l l J i(. r ')\ 2 to define the energy terms E s +E* — 
5 / | t rd r'| " 2Hx ( r ' r ') an d we se t the EXX total energy to 
e exx = j dr [_ij2 i ^( r )V 2 iP,(r)+n(r)V KS {r)]+E H + 
E*. 

Thus the total energy of a given system can be calcu- 
lated exactly from the KS potential, with the exception 
of one term: the correlation energy, defined here through 
E c = E — E EXX where E is the true groundstate en- 



ergy of the system. The correlation energy term essen- 
tially bundles the "difficult" physics of the true many- 
electron system into a single term, which is a highly non- 
local functional of the density and/or Kohn-Sham orbital 
wave functions, and must be approximated. An ab ini- 
tio way to evaluate correlation energies is to use time- 
dependent DFT [8] via the linear density-response func- 
tion, the fluctuation-dissipation theorem, and the adia- 
batic connection formula to form the "ACFD" functional. 
In recent years there has been a large increase in the use 
of ACFD functionals, particularly for the evaluation of 
vdW dispersion. The majority of these also make use of 
the direct random-phase approximation (dRPA) which 
we define later. A good discussion on, and summary 
of the ACFD-dRPA approach can be found in Ref. 9, 
although initial calculations on inhomogeneous systems 
were carried out more than a decade ago [10]. 

Theoretically exact applications of the ACFD involve 
the unknown dynamic exchange-correlation kernel / xc , 
a two point function defined as the the second func- 
tional derivative of the xc energy via / xc (r,r';i — t') = 
S 2 E XC /[6n(r,t)Sn(r' ,t')]. The concept of the xc ker- 
nel can also be extended to current-response theory 
where the tensor kernel F xc is known [11] to be a more 
'amenable' functional of the density. In practice / xc must 
be approximated, and the commonly employed dRPA in- 
volves setting / xc w 0. Perhaps surprisingly, the ACFD- 
dRPA functional has generally performed well for cal- 
culating energy differences, but not so well for abso- 
lute energies. Through the years various approxima- 
tions have been proposed for the / xc kernel, including 
the ALDA[12], energy-optimised kernel[13] and the Pe- 
tersilka, Gossman and Gross exchange kernel[14]. These 
have met with varying degrees of success in different sys- 
tems, but none has worked well in a wide range of sys- 
tems. More recently the exact exchange kernel / xc « / x 
has been evaluated[15-18] in the time-dependent EXX 
(tdEXX) approach leading, via the ACFD functional, 
to excellent results for correlation energies of atoms 
and molecules. However this kernel is very difficult 
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[0(N 5 )/0(N 6 ) in molecular basis function language] to 
calculate in practice, requiring inversion of the response 
or solutions of non-linear eigen-equations. Similarly, al- 
ternative approaches such as RPAx[19] and SOSEX[20] 
exist to improve on the ACFD-dRPA by including many- 
electron exchange but again these are numerically more 
difficult problems than the dRPA. 

The ISTLS formalism[21, 22], extending a total energy 
method for jellium[23] to general systems, was developed 
as a means of approximating the dynamic interactions in 
a sophisticated manner by making use of a self-consistent 
pair-correlation function. As shown in Rcf. 22 it is equiv- 
alent to self-consistently approximating F xc in an ACFD 
functional and it has so far enjoyed success in semi- 
homogeneous test systems [21, 24-26], most notably cor- 
rectly reproducing the difficult transition from a three- 
to a two-dimensional metal, something the dRPA fails to 
do. In some sense the ISTLS represents the 'next step' 
of ACFD-like approximation: introducing self-consistent 
physics to the dynamic tdDFT calculation in a rigorous 
manner through F xc , rather than deriving f xc or F xc from 
the groundstate calculation. 

In the original paper [21] on the method, the ISTLS 
functional was also tested on the helium atom where it 
performed very well, calculating the correlation energy to 
within O.lmHa. Advances in computing power and im- 
provements in numerical techniques have since allowed 
for wider testing. Here we discuss the implementation of 
the functional in spherical systems, and test it in a set 
of spherically symmetric neutral atoms and ions, includ- 
ing spin-polarised systems such as atomic sodium and 
lithium. 

The Kohn-Sham equations hipi = Ciipi can be used 
to generate the one-electron like orbitals of a system 
with a time-invariant KS potential V KS . In the ab- 
sence of a magnetic field but the presence of a small 
perurbation to V KS of form 5V(r,t) = <5U(r)e 1 "* we 
can write the change in density of the system as Sn — 
J dr'xo{r,r';uj)SV(r') = J dr'i/ (r,r';uj) ■ VSV(r') 
where xo is the bare (non-interacting) density-density 
response of the system, and v>q is the bare vector re- 
sponse. The change in current can be defined via 
Sj = iuj J dr'P (r, r'; uo)V5V(r') where P is the bare 
current-current response. Using tensor notation[27], it 
follows from these expressions that xo = — V' • u n and 
isq = — V • P. Each of these has an interacting equiva- 
lent eg. x\ which corresponds to the response a related 
system with electron-electron Coulomb interactions of 
strength A but with the groundstate density unchanged. 
When A = 1 these are equivalent to the response of the 
system to a change in the external potential. 

The ACFD correlation functional can be defined as 

E c = J dXj — J drdr'<S> x (r,r',ts) (1) 

with integrand [2 7] $ A (r,r',w) = [x\-Xo](r,r';w)v(\r'- 
r\) = [u x - v ](ry-,w) ■ V'v(\r' - r\) = [P A - 
P ](r,r';w) : V(|r' - r|). Here v(R) = 1/R is the 



Coulomb potential and V(i?) = —V <8> 'Vv(R) is its ten- 
sor equivalent. We can explicitly write the bare density- 
density and density-current reponses as 

Xo (r,r'; JS ) =2»^/ i ^(r)^(r , )G i (r,r , ) ) (2) 

i 

Mr, r'; is) -3 £ / 4 [t'(r)*(r')V'G,(r, r') 

i 

-Gi(ry)VVf(r^(0]/< (3) 

where Gi is short-hand for the bare one-electron Greens 
function G(r, r'; e,+is), a solution of [h — Sl]G(r, r'; 12) — 
S(r — r'). The current-current response Po has a similar 
expression. The interacting responses are defined via 

XA =Xo + Xo * (Aw + fx c ) * X\ (4) 
Pa =Po + Po*(AV + F^ c )*P A (5) 

where A* B = J dxA(r, x)B(x, r') and we take tensor 
products where appropriate. It is only in this relationship 
between the interacting and non-interacting case that the 
xc kernel is involved. 

The ISTLS scheme can be written[22] as a tensor F xc 
of form 

AV+F^l^OV^^V', (6) 

gx(r,r') ^n 2 x(r 7 r')/[n(r)n(r / )} (7) 

where n 2 \ is the interacting groundstate pair density 
at coupling strength A and n is the groundstate den- 
sity. Here we self-consistently calculate the dynamic 
interactions via the pair density n 2 \ calculated by the 
fluctuation-dissipation theorem 

n2\{i' 1 r') =n(r)n(r') — 5(r — r')n (r) 

-J ^XA(r,r'; JS ), (8) 

XA(r,r';*s)=(V®V') : P x (r,r'). (9) 

In practice we must iterate these equations: i) set g\ w go 
(ie. Hartree and exchange only) such that go{r,r') = 

1 - [n°(r)n°(r')]- 1 |E 4 / l V' 4 (r-)V'r(r-')| 2 , ") calculate P A 
via (5) and (6), iii) use P A to calculate a new g\ via (8) 
and iv) use the new g\ in ii) and repeat until convergence 
is reached. 

Making use of (6) and (9) we can transform (5) into 
Xa = Xo + Qx * Xa where Q x (r,r') = J dxv {r,x) ■ 
F\(x,r') and 

F x (ry)=g x (r,r')Vy^. (10) 

Thus it is possible to evaluate the ISTLS equations using 
only xo and Vq and not the full tensor current-current re- 
sponse P - This form of the equations is the original [21] 
approach to ISTLS calculations. It should be noted 
that the Petersilka-Gossman-Gross (PGG) kernel[14] 
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can be defined in a similar manner with Q\(r,r') = 
J dxu (r, x) ■ V x90 (x, r')^ ee J dr Xa (r, x)*f£$. 

In spherically symmetric atoms we can separate the 
orbitals as ipi(r) = ip n i m (r) = R n i(r)Yi m (r) and 
£i = where Y; TO is a spherical harmonic func- 
tion. The potential is V KS (r) = V KS (r) and the 
radial function satisfies hiR n i(r) = e n iR n i(r) where 
hi = -lir^drdrr - 1(1 + l)r- 2 } + V KS (r) and d r = 
d/dr. It follows from the properties of spherical 
harmonics and the definition of the Greens function 
that £ m C*r»Vw(r') = ^Pi{x) lnl (r,r') and 
G(r,r';Q) - Zi^Pi(x)G?(r,r') where x = r ■ r\ 
Pi (x) is a Legendre polynomial of order I and we use the 
short-hand J n l(r,r') — R n i(r)R n i(r'). Here Gp satisfies 
[hi - 0]Gp(r,r') = 5{r - r')/(rr'). It also follows from 
the symmetry of the system that 
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(11) 



X\(r,r';is) =^ — P L {x)xxL(r,r';is) 

L 

or i i 

u x (r, r>; is) = £ —±-P L { x )[v XL r' + «^r'J. (12) 

where r' ± = r — (r ■ r )r = r — xr . Thus the response 
equation is diagonal in L and xxl = Xol + Qxl *r Xxl 
where A * r B = J °° R 2 dRA(r, R)B(R, r'). 

Making use of the completeness of the polynomials 
Pi(x) we define the bare (A = 0) responses through 

X ol(t, r'; is) =2 ]T Kf vln m\V^ s (13) 

nil' 

ul L {r, r'; is) =± £ ^{ 7 „ ; [^3G^ ;+JS ] 

- [^7ni]^" i+4S } (14) 
ut L (r,r';is) =^ - ^n^'**. (15) 



mi' 

The Clebsch-Gordan-like coefficients K^, and fify, are 
defined as = &^2^1 J^dxPiP u P L and 

= ^^WtTTT !-i^DiP v P L where A = [d x Pi{x)\. We 
can similarly expand the vector kernel (10) of the ISTLS 
scheme as 

or i i 

F ( r ' r ') =Y,^ r PL{xW L {r,r')r + F^r ± ] (16) 

where F£ = Y.w K ii'9xi[d r vv] and Ft 
J2w Pvi9\ivv jr. We define g X i through (7) and (8) but 
with x\i{ r i r ') oru yi and use the Legendre expansion of 
the Coulomb potential l/\r — r'\ = J2 l vi(r, r') ( 2 ' + *) p '( x ) 
to define v t = min(r, r') 1 max(r, r')~(' +1 ) . Finally, 
using (12) and (16) we find[28] 

Qxl =v r 0L * r Fl + k[^ L * r Ft] - [k^ L ] *r [kFt] (17) 



We note that, with the exception of the self-consistency 
condition [defined via (8)], all terms are diagonal in 
s but couple together different I and involve convo- 
lutions over radial co-ordinate r. This allows us to 
evaluate Xxl{t, r'; is) from the sets {xoi( r , r '] is)}i and 
{voi(r, r'; is)}i provided the set {gxi}i is already known. 
Once Q\L{r,r';is) is calculated the system is diagonal 
in L and convolutions are only ever taken across r. In 
spin-polarised systems we must also introduce spin a =J \^ 
such that all radial coordinates are replaced by ra and 
convolutions include a sum over spin. 

To solve such a system numerically, we choose a grid 
of up to 512 radial points, and solve for the groundstate 
using the method of Krieger, Li and Iafrate[29] (KLI). 
The KLI approximation predicts E EXX quite accurately, 
and reproduces the correct — 1/r tail in atoms, a feature 
not present in LDA or GGA calculations. As such we 
feel it is an ideal starting point for these calculations. 

The grid {ri}, its weights {wi}, the radial orbital wave- 
functions R n i(ri), KS energies e„;, and the KS potential 
V KS (ri) are then stored for later use in the calculation 
of xo and vq- The Greens function can be solved quickly 
at arbitrary I and fi via a shooting method such that 



G?(r,r') 



1 \I(r)0(r') r<r' 



2rr'Wr 1 0(r)I(r') r>r' 



(18) 



where Wr = Id r O — Od r I and I(r) and 0(r) are "inner" 
or "outer" solutions of [hi —il]X(r) = with the bound- 
ary conditions I(r — > 0) cx r l and Oir — \ oo) = 0. Its 
radial derivative is then d r > Gp = Dp — Gp jr' where 



D?(r,r') 



1 



l{r)d r ,0(r') r<r' 
~2rr'Wr 1 0{r)d r >I{r') r > r' ' 



(19) 



where kf L = K^+\f L+1 + K^' h-i- 



L-l . 



We choose a set of abcissae and weights for s based 
on a Clenshaw-Curtis quadrature scheme, chosen for its 
accuracy in integrating Lorentz functions, such that con- 
vergence is reached using at most 50 points. We also 
exploit the fact that the system is diagonal in s to calcu- 
late and store response functions at a single s only and 
cumulatively evaluate integrals for the pair density and 
correlation energy. The method is also diagonal in A and 
we solve to high accuracy using A = ^ , | , 1 with appro- 
priate weights. We must also choose a cutoff in L which 
we set at L max = 6. In all tested cases the contribution 
to the energy from the L = 5 term is under 0.5%, with 
at least 97% of the energy accounted for by L < 3. 

Calculations are thus performed as follows: 1) for 
each I form the matrices XM{ri 7 rj;is) 7 ^(rj, r^; is) and 
vtl( r ii r j't * s ) and, at the first iteration, g\i w goi(ri,rj); 
2) take the stored response functions and pair densi- 
ties {g\i}i, then use quadrature to form QxL{fi,rj) via 
(17); 3) solve the matrix equation xxLij = XoLij + 
J^kQxLikWkXxLkj, repeating l)-3) for each L and each 
s 4) calculate new values for {g\i}i through a weighted 
mix of the existing data and the newly evaluated [via (8)] 
{gxi}f, 5) repeat from 1) until converged; 6) reset {g\i 
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and repeat from 1) for a new A. Typically it takes be- 
tween four and six iterations mixing 70% new and 30% 
old pair density to converge a correlation energy. It is 
worth noting that at each stage we impose symmetry 
under exchange of r and r' on each g\i- While formally 
this may differ slightly from the true ISTLS method, tests 
indicate that the correlation energy remains virtually un- 
changed, while convergence is improved. 



TABLE I. Correlation energies (in -mHa) for spherical atoms 
and ions. Includes the mean absolute error % (MAE%) for the 
neutral atoms (N), ions (I) and all atoms and ions together. 
He* is the extrapolation to Z — oo for a two-electron system. 



Atom 


RPA 


PGG 


ISTLS 


tdEXX" 2 


Exact 6 


He 


84.0 


44.9 


42.3 


44 


42.0 


Li 


113 


49 


41 


- 


45 


Be 


181 


104 


79 


102 


94 


N 


336 


145 


191 


- 


188 


Ne 


585 


331 


405 


389 


390 


Na 


612 


329 


413 




396 


Me 


672 


374 


458 


445 


438 


P 


833 


418 


563 




540 


Ar 


1071 


578 


744 


721 


722 


MAE% N 


76 


15 


5 






H" 


74.9 


43.6 


36.4 




42.0 


Li+ 


86.7 


45.3 


42.8 




43.1 


Be 2+ 


88.3 


45.6 


43.7 




44.3 


Ne 8+ 


91.1 


46.1 


45.4 




44.7 


Hg 78+ 


92.4 


46.3 


46.2 




46.5 


He* 


92.7 


46.4 


46.4 




46.9 


Be+ 


124 


51 


37 




47 


LP 


146 


84 


69 




73 


B+ 


207 


120 


86 




111 


Na+ 


582 


323 


404 




389 


Mg+ 


623 


331 


422 




400 


MAE% I 


94 


7 


7 






MAE% 


86 


11 


6 







a From Ref. 15, b From Refs 30-33 



In Table I we present correlation energies for a vari- 
ety of spherically symmetric systems. We compare the 
ISTLS energies with those from the dRPA and PGG cal- 
culated using the same code, with tdEXX energies from 
Ref. 15, and with 'exact' correlation energies from bench- 
mark methods [30-33]. We also include an extrapolation 
to the Z = oo case for the Helium isolectronic series (la- 
beled He*) by fitting E C (Z) vs. 1/Z for Z > 3. We have 
included only those atoms and ions that converged under 
the ISTLS self-consistency loop with a reasonable mix- 



ing coefficient and thus reasonable time. For benchmark- 
ing we compared our dRPA results with those of Jiang 
and Engel[31] and found agreement well within expected 
methodological bounds. 

In general the ISTLS does very well for correlation en- 
ergies, outperforming the dRPA in all tested systems, and 
the PGG in all but a few systems. In all the systems bar 
He where comparable tdEXX results were available[15] 
it outperforms the ISTLS, however this accuracy comes 
at much greater computational expense. ISTLS performs 
less well for ions than for atoms, with the greatest error in 
Bc + and B + . It is possible that, in these cases, the ISTLS 
iterations converge to an incorrect result, however test- 
ing this is difficult. For C 2+ the ISTLS method did not 
converge at all, most likely due to numerical instabilities 
in the high-density core region. It is worth noting that 
the ISTLS always pulls the PGG results back towards 
the true value, albeit overly so in some cases. While the 
PGG approximation performs slightly better than ISTLS 
for some of the smaller systems tested here, it is known 
to break down in bulk systems, particularly low density 
metals where it under-correlates[34]. This failure can be 
seen in the trend presented here, where the relative ab- 
solute PGG error increases with system size while ISTLS 
improves. By contrast the ISTLS performs consistently 
well for jcllium[23], metallic surface energies[24], across 
two- and three-dimensional metals[25], and here in the 
spherical atoms and ions. 

The numerical cost of the ISTLS functional scales with 
system size in a similar manner to standard ACFD- 
dRPA methods, but with a larger pre- factor and slightly 
larger memory requirements. In the best case sce- 
nario, the ISTLS can scale as (9(iV 4 ), while tdEXX and 
RPAx can scale as 0(N 5 ), a saving of one order. Our 
ISTLS calculations took between ten and twenty times 
as long as the ACFD-dRPA and used around five times 
the memory. The detailed method presented here may 
point the way to implementation in more general geome- 
tries involving expansions on Gaussian-type and Slater- 
type orbitals[35]. Implementation in existing plane-wave 
based bulk ACFD-dRPA codes should also be possible, 
albeit with non-trivial changes. 

Overall, we believe that the ISTLS is an excellent can- 
didate for a 'next step' functional, going beyond the 
dRPA. The tests on spherical systems further confirm 
its versatility, showing accurate results in systems with 
vastly different physics to those previously tested. With 
work on efficiencies and implementation it could, in fu- 
ture, provide viable benchmark calculations for electronic 
systems where existing high-level methods, such as the 
popular ACFD-dRPA, fail to achieve the desired level of 
accuracy and where wavefunction methods are too diffi- 
cult. 

The authors were supported by ARC Discovery Grant 
DP1096240. 
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